Impact on air-traffic -Geospatial

On major European routes

#load data on top European routes extracted from MongoDb
european_2020_routes <- read_csv("Data/top_european_routes_2020.csv")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_character(),
##   day = col_datetime(format = ""),
##   firstseen = col_datetime(format = ""),
##   lastseen = col_datetime(format = ""),
##   altitude_1 = col_double(),
##   altitude_2 = col_double(),
##   origin_lon = col_double(),
##   origin_lat = col_double(),
##   dest_lon = col_double(),
##   dest_lat = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
european_2020_routes %>% 
  group_by(origin) %>% 
  summarise(n())
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 22 x 2
##    origin `n()`
##    <chr>  <int>
##  1 EBBR    5321
##  2 EDDB    1937
##  3 EDDF   12688
##  4 EDDM    9814
##  5 EGCC    5051
##  6 EGLL   16750
##  7 EGPH    3206
##  8 EHAM   13277
##  9 EIDW    8146
## 10 ENGM    4710
## # … with 12 more rows
#filter out destinations that are not in Europe and convert dates using lubridate
european_2020_routes <-  european_2020_routes %>% 
  filter(origin!=destination) %>% 
  mutate(day= as_date(day, tz = NULL),
         route= gsub(" ", "", paste(origin,"-", destination)))%>% 
  filter(origin!="KJFK", destination!="KJFK") %>% 
  filter(day>"2020-01-31")

#load shape file for Europe
Europe <- ne_countries(scale = 'large', type = 'map_units', returnclass = 'sf', continent="Europe")

#trim map to get only western Europe
Europe <- sf::st_crop(Europe, xmin = -15, xmax = 40, ymin = 25, ymax = 62)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
#plot each flight in Europe between start of Feb and end of April
p <- ggplot() +
  geom_sf(data = Europe, size = 0.125) +
  geom_curve(
    data = european_2020_routes, 
    aes(x = origin_lon, y = origin_lat, xend = dest_lon , yend = dest_lat),
    curvature = 0.2, color= "#2c7fb8", size=0.1, alpha=0.4)+
  theme_void()+
  theme(
    plot.title = element_text(color="black", size=16, face="bold"),
    plot.subtitle = element_text(color="black", size=14, face="plain")
  )

p <- p+ transition_time(day) +
  labs(title = "Air traffic all but dissapeared in European airspace after COVID-19 struck",
       subtitle= "Daily flights between 21 major airports dropped sharply in late March \n\nDate: {format(frame_time, '%b %d, %Y')}"
  )

#Use gganimate to create an animated visualization
animate(p, nframes = 80, fps=7, height = 800, width =1200)

## On major and long global routes

#load data on global routes extracted from MongoDb
global_2020_routes <- read_csv("Data/top_global_routes_2020.csv")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_character(),
##   day = col_datetime(format = ""),
##   firstseen = col_datetime(format = ""),
##   lastseen = col_datetime(format = ""),
##   altitude_1 = col_double(),
##   altitude_2 = col_double(),
##   origin_lon = col_double(),
##   origin_lat = col_double(),
##   dest_lon = col_double(),
##   dest_lat = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
global_2020_routes <-  global_2020_routes %>% 
  filter(origin!=destination) %>% 
  mutate(day= as_date(day, tz = NULL),
         route= gsub(" ", "", paste(origin,"-", destination)))

#load shape file of the earth from rnaturalearthhires
world <- ne_countries(scale = "large", returnclass = "sf") %>%
  filter(name != "Antarctica") 

p <- ggplot() +
  geom_sf(data = world , size = 0.125) +
  geom_curve(
    data = global_2020_routes, 
    aes(x = origin_lon, y = origin_lat, xend = dest_lon , yend = dest_lat),
    curvature = 0.2, color= "#2c7fb8", size=0.1, alpha=0.4)+
  theme_void()+
  theme(
    plot.title = element_text(color="black", size=16, face="bold"),
    plot.subtitle = element_text(color="black", size=14, face="plain")
  )

p <- p+ transition_time(day) +
  labs(title = "COVID-19 hit long international routes very hard",
       subtitle= "Daily intercontinental flights between 31 major global airports \ndropped sharply in late March- early April  \n\nDate: {format(frame_time, '%b %d, %Y')}"
  )

#Use gganimate to create an animated visualization
animate(p, nframes = 80, fps=7, height = 800, width =1200)

Virus vs the 2010 Volcano

#load data on 7D moving averages of flights for 2010,2019, 2020
air_crises <- read_csv("Data/COVID_7D_Time series.csv") %>% 
  #clean column names
  janitor::clean_names() %>% 
  rename(cases=daily_new_covid_19_cases_7_day_moving_avg)
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   Date = col_character(),
##   `2010` = col_double(),
##   `2019` = col_double(),
##   `2020` = col_double(),
##   `Daily new COVID-19 cases (7 day moving avg.)` = col_double()
## )
#skim to view the data
skimr::skim(air_crises)
Data summary
Name air_crises
Number of rows 366
Number of columns 5
_______________________
Column type frequency:
character 1
numeric 4
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
date 0 1 11 12 0 366 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
x2010 1 1.00 26180.98 3349.78 12372 23734.0 26796.0 29183.0 30221.0 ▁▁▂▅▇
x2019 1 1.00 30122.22 3796.30 22639 26469.0 30265.0 34204.0 35235.0 ▂▇▂▅▇
x2020 61 0.83 14263.26 7715.80 3051 6234.0 14566.0 20506.0 26352.0 ▆▁▇▁▆
cases 61 0.83 27046.80 40811.00 0 7196.4 12511.4 31298.4 233806.6 ▇▁▁▁▁ 
air_crises <- air_crises %>% 
  mutate(
    date= as.Date(date, format= "%b %d,%Y")
  ) %>% 
  select(-cases) %>% 
  #pivot longer for to facilitate visualisation in ggplot
  pivot_longer(cols= "x2010":"x2020", names_to = "line", values_to = "value") %>% 
  filter(date<="2020-09-30")


ggplot(air_crises, aes(x=date, y=value, color=line))+
  geom_line()+
  theme_minimal()+
  theme(
    axis.title.y = element_blank(),
    axis.title.x =element_blank(),
    axis.text.x = element_text(color="black", size=12, face="plain"),
    axis.text.y= element_text(color="black", size=12, face="plain"),
    plot.title = element_text(color="black", size=20, face="bold"),
    plot.subtitle = element_text(color="black", size=16, face="plain"),
    legend.text = element_text(size = 12), 
    legend.title = element_blank()
  )+
  labs(
    title="Virus much worse than Volcano!",
    subtitle= "7 day moving average of number of daily flights in EuroControl member states\n\n",
    color="Year"
  )+
  #define color and labels manually
  scale_color_manual(values= c("#3288bd","#1a9850" ,"#762a83"), labels = c("2010", "2019", "2020"))+
  #adding a text annotation to draw attention to the 2010 crises
  annotate("text", x= as.Date("2020-05-20"), y=20000, label= "Air travel disruption \nafter 2010 eruptions of \nEyjafjallajökull ") +
  #cleaning up 0's on the y axis
  scale_y_continuous(breaks=c(10000,20000,30000), labels= c("10k","20k","30k"))

Impact on vertical flight efficiency and additional taxi-out time

#load data on flight efficiency during descent
flight_efficiency_desc<- read_csv("Data/COVID_CDO_Page 1_Time series.csv") %>% 
  clean_names() %>% 
  mutate(
    date_2= as.Date(date_2, format= "%b %d,%Y"),
    year=as.factor(year)
  ) 
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   DATE_2 = col_character(),
##   YEAR = col_double(),
##   `Avg. time in level flight (seconds)` = col_double()
## )
#load data on flight efficiency during ascent
flight_efficiency_asc<- read_csv("Data/COVID_CCO_Page 1_Time series.csv") %>% 
  clean_names() %>% 
  mutate(
    date_2= as.Date(date_2, format= "%b %d,%Y"),
    year=as.factor(year)
  ) 
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   DATE_2 = col_character(),
##   YEAR = col_double(),
##   `Avg. time in level flight (seconds)` = col_double()
## )
#calculate 99% confidence intervals for mean time level during CDO
flight_efficiency_plot <- flight_efficiency_desc %>% 
  group_by(year) %>% 
  summarise(avg_time_level_desc= mean(avg_time_in_level_flight_seconds,na.rm=TRUE),
            sd_time_level_desc= sd(avg_time_in_level_flight_seconds, na.rm=TRUE),
            count = n(),
            se_time_level_desc = sd_time_level_desc/sqrt(count),
            #compute t-critical value
            t_critical = qt(0.99, count-1), 
            margin_of_error = t_critical * se_time_level_desc,
            #calculate upper and lower bound of the confidence interval
            lower_bound = round((avg_time_level_desc - margin_of_error),digits=2),
            upper_bound = round((avg_time_level_desc+ margin_of_error), digits=2)
  ) 
## `summarise()` ungrouping output (override with `.groups` argument)
#plot the 99% confidence intervals for mean time level during CDO
p1 <- ggplot(flight_efficiency_plot, aes(x=year, y=avg_time_level_desc, colour=year)) +
  geom_point(size=3) +
  #create the error bar using the confidence intervals calculated in previous block
  geom_errorbar(width=.5, aes(ymin=lower_bound, ymax=upper_bound),  size= 1) + 
  labs( y= "Mean time in level flight (seconds)", 
        title="Vertical flight efficiency is significantly higher in 2020",
        subtitle= "99% confidence intervals of mean time in level flight across EuroControl member states \n\nDuring continuous descent operations (Landing)"
  ) + 
  theme_minimal()+
  coord_flip()+
  theme(
    #clean up text sizes of various elements in the graph
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 13),
    axis.text.x = element_text(color="black", size=12, face="plain"),
    axis.text.y= element_text(color="black", size=12, face="plain"),
    plot.title = element_text(color="black", size=16, face="bold"),
    plot.subtitle = element_text(color="black", size=14, face="plain"),
    legend.text = element_text(size = 12), 
    legend.title = element_text(size = 13),
    legend.position = "none"
  )+
  scale_color_manual(values= c("#636363","#2c7fb8"))

#calculate 99% confidence intervals for mean time level during CCO
flight_efficiency_plot <- flight_efficiency_asc %>% 
  group_by(year) %>% 
  summarise(avg_time_level_asc= mean(avg_time_in_level_flight_seconds,na.rm=TRUE),
            sd_time_level_asc= sd(avg_time_in_level_flight_seconds, na.rm=TRUE),
            count = n(),
            se_time_level_asc = sd_time_level_asc/sqrt(count),
            #compute t-critical value and margin of error
            t_critical = qt(0.99, count-1), 
            margin_of_error = t_critical * se_time_level_asc,
            lower_bound = round((avg_time_level_asc - t_critical * se_time_level_asc),digits=2),
            upper_bound = round((avg_time_level_asc + t_critical * se_time_level_asc), digits=2)
  ) 
## `summarise()` ungrouping output (override with `.groups` argument)
#plot the 99% confidence intervals for mean time level during CCO
p2 <- ggplot(flight_efficiency_plot, aes(x=year, y=avg_time_level_asc, colour=year)) +
  geom_point(size=3) +
  #create the error bar using the confidence intervals calculated in previous block
  geom_errorbar(width=.5, aes(ymin=lower_bound, ymax=upper_bound),  size= 1) + 
  labs( y= "Mean time in level flight (seconds)",
        title="",
        subtitle= "\n\nDuring continuous climb operations (Takeoff)",
        caption="Source: Eurocontrol Network"
  )+
  theme_minimal()+
  coord_flip()+
  theme(
    #clean up text sizes of various elements in the graph
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 13),
    axis.text.x = element_text(color="black", size=12, face="plain"),
    axis.text.y= element_text(color="black", size=12, face="plain"),
    plot.title = element_text(color="black", size=16, face="bold"),
    plot.subtitle = element_text(color="black", size=14, face="plain"),
    legend.text = element_text(size = 12), 
    legend.title = element_text(size = 13),
    legend.position = "none"
  )+
  scale_color_manual(values= c("#636363","#2c7fb8"))

#place plots one on top of the other using patchwork
p1/p2

#Create time-series plots of mean time level during CCO and CDO
# To read up on Continuous climb and descent operations visit
#https://www.eurocontrol.int/concept/continuous-climb-and-descent-operations

#Tine series plot for CDO
ggplot(flight_efficiency_desc, aes(x=date_2, y=avg_time_in_level_flight_seconds, color=as.factor(year)),size=2)+
  geom_line()+
  theme_minimal()+
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(color="black", size=12, face="plain"),
    axis.text.y= element_text(color="black", size=12, face="plain"),
    plot.title = element_text(color="black", size=16, face="bold"),
    plot.subtitle = element_text(color="black", size=14, face="plain"),
    legend.text = element_text(size = 12), 
    legend.title = element_blank()
  )+
  labs(
    title="COVID-19 improved vertical flight efficieny during Continuous Descent Operations (CDO)",
    subtitle= "Average time in level flight (in seconds) for aircrafts in EuroControl member states",
    caption="Source: Eurocontrol Network"
  )+
  scale_color_manual(values= c("#636363","#2c7fb8"), labels = c("2019", "2020"))

#Tine series plot for CCO
ggplot(flight_efficiency_asc, aes(x=date_2, y=avg_time_in_level_flight_seconds, color=as.factor(year)), size=2)+
  geom_line()+
  theme_minimal()+
  theme(
    axis.title.y = element_blank(),
    axis.title.x = element_blank(),
    axis.text.x = element_text(color="black", size=12, face="plain"),
    axis.text.y= element_text(color="black", size=12, face="plain"),
    plot.title = element_text(color="black", size=16, face="bold"),
    plot.subtitle = element_text(color="black", size=14, face="plain"),
    legend.text = element_text(size = 12), 
    legend.title = element_blank()
  )+
  labs(
    title="COVID-19 improved vertical flight efficieny during Continuous Climb Operations (CCO)",
    subtitle= "Average time in level flight (in seconds) for aircrafts in EuroControl member states",
    caption="Source: Eurocontrol Network"
  )+
  scale_color_manual(values= c("#636363","#2c7fb8"), labels = c("2019", "2020"))

#load data on additional taxi-out time
taxi_out_time<- read_excel("Data/Taxi Out Time.xlsx") %>% 
  clean_names() %>% 
  mutate(year=as.factor(year),
         #calculate average additional taxi-out time
         avg_unimpeded_taxi= time_txo_unimp_2/ flt_txo_unimp_2,
         avg_additional_taxi= time_txo_add_2/ flt_txo_unimp_2) %>% 
  filter(month_num>=3, month_num<=9, year %in% c(2018,2019,2020))

#compute 99% confidence intervals
taxi_out_time_plot <- taxi_out_time %>% 
  group_by(year) %>% 
  summarise(avg_add_taxi= mean(avg_additional_taxi,na.rm=TRUE),
            sd_add_taxi= sd(avg_additional_taxi, na.rm=TRUE),
            count = n(),
            se_add_taxi = sd_add_taxi/sqrt(count),
            t_critical = qt(0.999, count-1), 
            margin_of_error = t_critical * se_add_taxi,
            #calculate lower and upper bounds using the formula
            lower_bound = round((avg_add_taxi - t_critical * se_add_taxi),digits=2),
            upper_bound = round((avg_add_taxi + t_critical * se_add_taxi), digits=2)
  ) %>% 
  mutate(is_2020= as.factor(ifelse(year==2020,1,0)))
## `summarise()` ungrouping output (override with `.groups` argument)
#display 99% confidence intervals visually to show difference between taxi-out times in 2020 vs 2019
ggplot(taxi_out_time_plot, aes(x=year, y=avg_add_taxi, colour=is_2020)) +
  geom_point(size=3) +
  geom_errorbar(width=.5, aes(ymin=lower_bound, ymax=upper_bound),  size= 1) + 
  #add titles and subtitles using "\n" to space and align text
  labs( y= "Mean additional taxi-out time (minutes)", 
        title="Flights took off significantly faster in 2020",
        subtitle= "99% confidence intervals of mean additional taxi-out time across \nEuroControl member states (Between March and September)",
        caption= "Source: Eurocontrol Network"
  ) + 
  theme_minimal()+
  coord_flip()+
  theme(
    plot.title = element_text(color="black", size=16, face="bold"),
    plot.subtitle = element_text(color="black", size=14, face="plain"),
    legend.text = element_text(size = 12), 
    legend.title = element_text(size = 13),
    legend.position = "none",
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = 13),
    axis.text.x = element_text(color="black", size=12, face="plain"),
    axis.text.y= element_text(color="black", size=12, face="plain")
  )+
  scale_color_manual(values= c("#636363","#2c7fb8"))

taxi_out_time_apt <- taxi_out_time %>% 
  filter(month_num>=6, month_num<=9, year %in% c(2019,2020)) %>% 
  pivot_wider(names_from = "year", values_from= "avg_additional_taxi") %>% 
  mutate(airport_w_code= paste(apt_name, "-", apt_icao)) %>% 
  group_by(airport_w_code, month_num) %>% 
  rename(pre_covid_t="2019",
         covid_t="2020") %>% 
  summarise( pre_covid_t = max(pre_covid_t, na.rm=TRUE),
             covid_t = max(covid_t, na.rm=TRUE)
  ) %>% 
  filter(is.finite(covid_t), is.finite(pre_covid_t)) %>% 
  group_by(airport_w_code) %>% 
  summarise(mean_covid_t= mean(covid_t),
            mean_pre_covid_t= mean(pre_covid_t),
            count=n(),
            percentage_change= (mean_covid_t- mean_pre_covid_t)/mean_pre_covid_t
  ) %>% 
  filter(count==4) %>% 
  arrange(percentage_change) %>% 
  head(20)


ggplot(taxi_out_time_apt, aes(x=fct_reorder(airport_w_code, percentage_change, .desc = TRUE), y=percentage_change))+
  geom_col(fill= "#2c7fb8")+
  coord_flip()+
  labs(
    title="Top 20 airports that saw the largest decrease in additional taxi-out time",
    subtitle= "Percentage change in mean additional taxi-out time between 2020 and 2019 (June to September) \n",
    caption= "Source: Eurocontrol Network"
  ) + 
  theme_minimal()+
  theme(
    axis.title.y = element_blank(),
    axis.title.x = element_blank(),
    plot.title = element_text(color="black", size=16, face="bold"),
    plot.subtitle = element_text(color="black", size=14, face="plain"),
    axis.text.x = element_text(color="black", size=12, face="plain"),
    axis.text.y= element_text(color="black", size=12, face="plain")
  )+
  scale_y_continuous(labels = scales::percent)